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Abstract — This paper is devoted to adaptive long autoregres- 
sive spectral analysis when (i) very few data are available, 
(ii) information does exist beforehand concerning the spectral 
smoothness and time continuity of the analyzed signals. The 
contribution is founded on two papers by Kitagawa & Ger- 
sch [1,2]. The first one deals with spectral smoothness, in the 
regularization framework, while the second one is devoted to 
time continuity, in the Kalman formalism. The present paper 
proposes an original synthesis of the two contributions: a new 
regularized criterion is introduced that takes both information 
into account. The criterion is efficiently optimized by a Kalman 
smoother. One of the major features of the method is that it is 
entirely unsupervised: the problem of automatically adjusting 
the hyperparameters that balance data-based vs prior-based 
information is solved by maximum likelihood. The improvement 
is quantified in the filed of meteorological radar. 

Index Terms — Adaptive spectral analysis, long autoregressive 
model, spectral smoothness, time continuity, regularization, hy- 
perparameter estimation, maximum likelihood, meteorological 
Doppler radar. 

I. Introduction 

ADAPTIVE SPECTRAL ANALYSIS and time-frequency 
analysis are of major importance in fields as widely var- 
ied as speech processing [3], acoustical attenuation measure- 
ments [4,5], ultrasonic Doppler velocimetry [6], or Doppler 
radars [7-11]. Reference [12] gives a synthesis of the vari- 
ous methods for these problems, and provides a number of 
bibliographical introductions. 

The present paper focuses on short-time analysis: typically, 
for analysis of pulsed Doppler signals only 8 or 16 samples 
are available to estimate one spectrum, with possibly various 
shapes (multimodal or not, of large spectral width or not, 
mixed clutter, etc. . . ). Under such circumstances, the construc- 
tion of the sought spectra becomes extremely tricky on the 
sole basis of the samples. As a point of reference, let us recall 
that several hundred samples are usually needed to compute 
an averaged periodogram with a fair bias-variance compro- 
mise [13,14]. So, parametric methods have generally been 
preferred, among which autoregressive (AR) play a central 
role. The AR coefficients estimation is usually tackled in the 
Least Squares (LS) framework [15, 16]. These methods often 
provide a solution at points where non-parametric methods 
are useless; but when the number of data is very low, these 
techniques become, in their turn, useless, especially if various 
spectral shapes are expected due to model order limitations. 
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In order to construct a reliable image, structural information 
about the sought spectrum sequence must be accounted for 
Our investigation is therefore restricted to the cases in which 
two kinds of information are foreknown: spectral smoothness 
and time continuity. This a priori information is the foundation 
of the proposed construction. 

In the framework of stationary AR analysis, Kitagawa & 
Gersch proposed a method integrating the idea of spectral 
smoothness [1] by which a high-order AR model can be 
robustly estimated, thereby getting around the difficult prob- 
lem of order selection, and providing capability to estimate 
various spectral shapes. For the non-stationary case, and aside 
from [1], the same authors introduced in [2] a Markovian 
model for the regressor sequence in the Kalman formalism, in 
order to reflect time continuity. The present paper reviews [1] 
and [2] and makes an original synthesis suited to the special 
configuration of Doppler signals. A new Regularized Least 
Squares (RegLS) criterion simultaneously includes the spectral 
and time information and is optimized by a Kalman Smoother 
(KS). 

One of the major features of the method is that it is entirely 
unsupervised: the adjustment of parameters that weight the 
relative contributions of the observation versus the a priori 
knowledge is automatically set by maximum likelihood (ML). 

A comparative study is proposed in the context of 
pulsed Doppler radars. Special attention is payed to atmo- 
spheric and/or meteorological context imaging or identifica- 
tion: ground clutter, rain clutter, sea echos, etc. . . Adaptive 
spectral estimation of mixed clutter is achieved by means of 
several usual AR methods and the proposed one. The latter 
achieves qualitative and quantitative improvements w.rt. usual 
methods. 

The paper is organized as follows. Section mainly intro- 
duces notations and problem statement. Section Hill focuses on 
usual LS methods and usual adaptive extensions. The proposed 
method is presented in Section |IV] and Section |V] deals with 
the KS. The problem of automatic parameter estimation is 
addressed in Section [Vl] Simulation results are presented in 
Section IVIII Finally, conclusions and perspectives for future 
works are presented in Section IVIIII 

II. Problem statement 

The problem is that of processing pulsed Doppler signals 
from electronic scanning radars or ultrasound velocimeter 
The reader may consult [6, 7] for a technological review. The 
pulsed Doppler systems are such that the observed signals do 
not occur in the usual form of time-frequency problems. So, 
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Fig. 1. Simulated observations over 110 range bins with 8 samples per bin (corresponding to 8 Doppler pulses). The Ihs figure shows the true spectra 
sequence. The narrow zero-mean spectra characterizes ground clutter (bin 15 to 57). Rain clutter induces more or less broad, single-mode spectra (bin 35 to 
75). Lastly, sea echos resulting from wave phenomena exhibit two maxima (bin 56 to 95). The middle figure shows the real part and imaginary part of the 
data and the rhs one shows the associated peiiodograms. 



neither the usual time-frequency methods nor the one proposed 
by Kitagawa & Gersch can be directly applied, and part of the 
presented work consists in constructing an appropriate method 
for the encountered configuration. 

The measurements are available as a set of complex signals 
y — [yi, . . . , Um], depth-wise juxtaposed in M range bins. It 
is assumed that each = [Umi, ■ ■ ■ ,?/mAf]* is a iV sample 
vector extracted from a zero-mean stationary process. Fig. [T] 
gives a Gaussian simulated example over M = 110 bins for 
which N — 8 samples are observed per bin. The successive 
regressors are denoted a,„ = [a„ip], where m indicates the 
considered bin (to e N^^ = {1,2,..., M}) and p the order 
of the autoregression coefficient (p G Np). Let us note 
A = [ai, . . . , ajv/] S <C^^^ the collection of the whole set of 
coefficients. Let us also introduce r„i and for signal and 
prediction error powers. The remainder of the paper is devoted 
to estimation of these quantities. The next section deals with 
the usual LS methods and their adaptive extension, and shows 
their inadequacy for the problem at stake. 

III. Review of classical methods 

A. Stationary spectral analysis 

This subsection is devoted to spectral analysis applied to 
a single bin to. Assuming a Gaussian distribution for the ob- 
served signal, the likelihood of the AR coefficients / {y,n\o.m) 
shows a special form [17, p. 82], but its maximization raises 
a difficult problem. A few authors [18, 19] have undertaken to 
solve it; but, firstly, the available algorithms cannot guarantee 
global maximization, and secondly, they are not computation- 
ally efficient for the applications under the scope of the paper 
To remedy these disadvantages, the following approximation 
of the likelihood function is usually accepted [16, p. 185]: 

/ (ymlflm) = (ttC)-^ exp (-Q^'(a„)/r^J , (1) 



involving the norm of the prediction error vector 



(2) 



i.e., a quadratic form w.r.t. the a„j namely the LS criterion. The 
Um and Ym are the vector and matrix designed according to 



some chosen windowing assumption [15, p. 217], [20, eq. (2)]. 
There are four possible forms: non-windowed (covariance 
method), pre-windowed, post-windowed, double-windowed 
i.e., pre- and post- windowed (autocorrelation method). Let us 
note L the size of y,„: L = N-P, L = N 01 L = N + P, 
according to the chosen form. This choice is of importance 
since it strongly influences spectral resolution for short time 
analysis [15, p. 225]. 

Whatever the chosen form, the maximization of ([T]i comes 
down to the minimization of (|2]l and yields: 



argminQJ;^(a™) = {Y^Yrr, 



-'Yl 



Vr, 



(3) 



As a prerequisite, the problem of choosing the model order 
P must be tackled: P has to be high enough to describe various 
PSD, and low enough to avoid spurious peaks, i.e., to ensure 
spectral smoothness. This compromise can usually be set by 
means of criteria such as FPE [21], AIC [22], CAT [23], or 
MDL [24], but, in the situation of prime interest here, they 
fail because the available amount of data is too small [25]. 
Actually, there exists no satisfying compromise in term of 
model order, since too few data are available to estimate DSPs 
with possibly complex structures. 

B. Adaptive spectral analysis 

For the "multi range bin" analysis, the first idea consists 
in processing each bin independently: according to the LS 
approach, it amounts to minimize a global LS criterion: 



(4) 



m— 1 



However, the resulting spectra hold unrealistic variations in 
the spatial direction (see Fig. |4]i. In order to remedy this 
problem, the Adaptive Least Squares (ALS) approach accounts 
for spatial continuity by processing the data from several 
bins, possibly in weighted form, to estimate each a,n- A first 
approach uses a series of LS criteria including the data in a 
spatial window of length W. A widely used alternative is the 
exponential decay memory which uses geometrically weighted 
LS criteria, with parameter A G [0, 1]. The latter is more 
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popular because it is simpler: A is merely incorporated into 
a standard recursive LS algorithm [15, p. 266]. In both cases, 
the degree of adaptivity, i.e., the spatial continuity is modulated 
by W or A. 

C. Conclusion 

Whatever the variant, the main disadvantage of these ap- 
proaches has to do with the parameter settings. 

- From the spectral standpoint, smoothness is introduced in 
a roundabout fashion, via the model order (adjusted by 
P) and the compromise no longer exists when the amount 
of data is reduced. 

- From the spatial standpoint, continuity is also indirectly 
introduced (and tuned by W or A) and no automatic 
method for adjusting this parameters is available. 

These hmitations are unavoidable in the simple LS formal- 
ism, and to alleviate this problem we resort to the regulariza- 
tion theory. In this framework, the proposed approach 

• includes the spectral smoothness and spatial continuity in 
the estimation criterion itself; 

• allows long-AR model to be robustly estimated, and then 
various spectra to be identified; 

• provides automatic parameter setting, i.e., an entirely 
unsupervised method. 

IV. Long AR - spatial continuity - spectral 

SMOOTHNESS 

A. Spatial continuity model 

The first idea consists in building a spectral distance. 
Following [2], starting with the PSD in bin m 

P 



(5) 



the proposed spectral distance between Sm and Sm' is founded 
on the fc-th Sobolev distance between A„-, and Am'- 



Dkijn, m') cx 



[A^{v) ~ A^.iv)] 



dv . 



Calculations similar to those of [2] yield a quadratic form: 

Dk{m,m') ^ (a,„ - a.^'^ Ak{a„i - am') , (6) 
where Afe = diag[l^'^, . . . , P^fcj jjjg spectral matrix. 

B. Spectral smoothness model 

The spectral smoothness measure proposed by Kitagawa & 
Gersch in [2] (see also [26]), is easily deduced from ^ as 
the distance to a constant DSP 

Dk{m) oc al^^kam ■ (7) 

According to [1,2], k £ 7L+, but Afe as well as Q and ^ 
can be extended to A; G IR+. 

Remark 1 — Strictly speaking, Dk{m,m') and Dk{m) are 
not spectral distances nor spectral smoothness measures since 
they are not functions of the PSD itself. However, they 
are quadratic and this has two advantages: it considerably 
simplifies regressor calculations (see Section as well as 
regularization parameter estimation ( see Section [Wl ). 



C. Double smoothness 

Starting with the spectral smoothness (|7]l and the spatial 
distance (|6]l, a new quadratic penalization is introduced: 



M 



M-1 



Q~(A) = - J"Dk{m) + — VDfe(m,m + l). (8) 

m—l m—l 

It integrates both spectral smoothness and spatial continuity 
respectively tuned hy Xg — l/r^ and Ad = l/r^. 

Remark 2 — The penalization dSj has a Bayesian interpre- 
tation [27] as a Gaussian prior for the sought regressors: 

/(A) cxGxp[-Q~(A)] , (9) 

useful for hyperparameter estimation, in Section 

D. Regularized least squares 

From the LS criteria (|4|i and the penahzation term the 
proposed RegLS criterion reads: 

g"-(A) - Q^%A)+Q^{A) (10) 
M ^ 

m=l 

M 



m—l 
^ M-1 

Tj ^ ^ 



involving three terms which respectively measure fidelity 
to the data, spectral smoothness and spatial regularity. The 
regularized solution is defined as the minimizer of ( fTOb : 



Anas = argmin(5''°'=(A) 

A 



(11) 



Remark 3 — The regularized criterion UOt has a clear 
Bayesian interpretation [27]: likelihood (O and prior ^ can 
be fused thanks to the Bayes rule, into a Gaussian posterior 
law for the sought regressors: 



f{A\y)^cxp[-Q^^-{A)] , 
So, the solution (|77} is also the MAP estimate. 



(12) 



E. Optimization stage 

Several options are available to compute ( fTTT i. Since 
Q^°^{A) is quadratic, -Aros is the solution of a AlP x AIP 
linear system. Moreover, since the involved matrix is sparse, 
direct inversion should be tractable but not recommendable 
here (M = 110, P = 7). Another approach may be found 
in gradient or relaxation methods [28] since (5"°*(A.) is 
differentiable and convex. But, given the depth-wise structure, 
another algorithm is preferred: Kalman Smoothing (KS). So, 
we resort to the initial viewpoint of Kitagawa & Gersch 
in [2]. However, it is noticeable that [2] does not mention 
the minimized criterion, where as our KS is designed to 
minimize (fTOl i. 
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V. Kalman smoothing 

A. State-space form 

- The successive prediction vectors are related by a 
first-order state equation: 



l-m+l 



(13) 



in which each £,„ is a complex, zero-mean, circular, 
vector with covariance matrix P^^^ — r^^^A^^ and the 
£,„-sequence, is depth-wise white. 

- The full state model also brings in the initial mean 
and covariance: the null vector and P° = r^A^T^, 
respectively. 

- The observation equation is the recurrence equation for 
the AR model in each bin, written in compact form as 



(14) 



i.e., a generalized version of the one proposed in [2], 
adapted to depth-wise vectorial data. Each e,„ is a com- 
plex, zero-mean, circular, vector with covariance r'^-^Jv, 
the e„j-sequence, is also depth-wise white. 

Remark 4 — [2] accounts for spatial continuity by means of 
a special case of Eq. ( I73I )." flm+i = ctm + £m- The latter has 
two drawbacks, though. Firstly, it is introduced apart from the 
idea of spectral smoothness. Secondly, from a Bayesian point 
of view, this equation is interpreted as Brownian process with 
an increasing variance, which may cause drifts to appear in 
the estimated spectra. On the contrary, the new coefficients am 
can be chosen in order to ensure stationarity of the model ( I73l > 
or to minimize the homogeneous criterion ( 1701 ). 

B. Equivalence between parameter settings 

1) Homogeneous criterion: This section establishes the 
formal link between the parameters of the KS (r" and rf^) 
and those of the regularized criterion ( fTOl l (rj and rg). [29, 
p.150-158] states that the KS associated to (O-lO mini- 
mizes: Q^^{A) 

M 



m— 1 
M-1 



re 



(15) 



Partial expansions yield identification of ( fTol ) and ( fTSI ) through 
the following count-down recursion. 
© Initialization (m = M — 1): 

ctM-i = (1 + , and r|,j_i = rd^M-i ■ 
(D Count-down recursion {m = M — 2, . . . , 1): 

a,„ = (2 + p - a„i+iy^ , and r";,-^ = rda,„ . 
@ The last step yields the initial power: 



with p = Td/r's e M*^. These equations allow to precompute 
the coefficients of the KS in order to minimize ( fTol l. 

2) Limit model: This section is devoted to the asymptotic 
behavior of the ckm-sequence. For the sake of notational 
simplicity, the sequence is rewritten in a count-up form: 



TO = 1 

TO e N* 



5i = (l + p)"' 



(16) 



It is clear that ai g]0,1[ since p E WC^. Let us introduce 
f{u) = (2 + p - uy^. It is straightforward that /(]0, 1[) C 
]0, 1[, so the entire Q;,„-sequence remains in ]0, 1[. Moreover, if 
it exists, the limit G [0, 1] necessarily fulfills f{a^) = a^. 
Elementary algebra yields: 



- Ve^ - 4) /2 



(17) 



with 9 = 2 + p = 2 + r^/r^. Finally, one can effortless see 
that V-u, V e]0, 1[ we have - /(t;)| < (1 + p)^^\u - v\, 

i.e., / is a Lipschitz function with ratio in ]0, 1[. Hence, the 
sequence effectively converges towards a^. It is also easy to 
see that the sequence is monotonous: increasing if ai < 
and decreasing otherwise. In the present case, comparison of 
ai in (fT6] l and a.^ in (fTTI l shows that the cim-sequence is 
decreasing (in the count-up form), hence, am is increasing. 

Finally, since rf„ = r^am, the corresponding limit state 
power is given by: 

= raa^ . (18) 

3) Associated stationary criterion: This section is devoted 
to the stationary limit model: the special case of Eq. ( fT3] l. 
with a„i — a^ and r^j = r^, i.e., a stationary first-order AR 
model for the a,„-sequence. The initial power is denoted r'^ 
for notational coherence, even if it is not defined as a limit. It 
is actually defined according to and a^ in order to ensure 



stationarity for the first-order AR model: r'^ 



7(1- 



a: 



Replacement of am , ^fn : "^oo : ''L i ''JL in Eq. (flST l yields 
the criterion minimized by the stationary KS: 



1 

Q (^) — ^ ^ ~iy^^ ^m^m)^ iym ^mO'n 
m=l ^™ 

jl-a^) 



M-l 



^ {O'm — am+l)^Afe(aTO — flm+l) 



a^{\ — a^j 



(a|Afeai + aj^jAfeO 



where superscript "S" stands for stationary. Since we have: 
Td = rl^/a^ from Eq. ( fTSl l and rg = 7'^/(l — ctoo)^ from 
Eq. ([TtI i. one can effortless see that: 

Q^A) - Q--(A) + "°°^^^7"°°^ alA,.ai + a\,AkaM) ■ 

So, the stationary criterion Q^{A) and the initial homogeneous 
one Q^°'^{A) are equal apart from the edge effects, i.e., two 
terms regarding the first and last regressors. As a consequence. 
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the minimizer of Q^'^'-IA) and Q^{A) are practically equiv- 
alent and the latter is preferred since it does not require 
precomputation of the and r^. 

C. Kalman smoother equations 
• Initialization (m = 1) 



ai\i = 

All = CA- 

Filtering phase (for m — 2, . . . , M) 
- Prediction step 



a 



m\va~l 



(19) 
(20) 



(21) 



Prn\m-i = aLPm-i\^-i+rl^k^ (22) 
- Correction step 





= P 1 

m m- 




(23) 








(24) 




= Vm - 




(25) 




^m rn- 


— 1 ^~ ^mP-ra ^ra 


(26) 




— P 1 


_1 — KjnRm K^-^ 


(27) 



• Smoothing count-down phase (for m = M — 1, . . . , 1) 

0„ = a.F™|™P„-ii|„ (28) 

0-m\M — Q-mlm + Qm {o-m+l\M ^ (^m+l\m) (29) 
,|M — Pm\m + Qm (^'m+l|M ~ Pm+l\rri) Qt?0) 



D. Fast algorithm 

Fast algorithms used to take a primordial position in the past 
decades, especially for real-time computations. More specif- 
ically, for adaptive spectral analysis of ultrasound Doppler 
signal, the Marasca algorithm [27] has been used in a real- 
time high-resolution velocimeter prototype. But, it has two 
drawbacks, resulting in a rigid spectral and spatial continuity 
tuning. On the one hand, it proceeds by blocks and incorpo- 
rates spatial continuity by using the regressor of the current 
block as a prior mean for the next one; on the other hand, the 
fast version is developed only for the zero-order smoothness 
(fc = 0). 

To our knowledge, no fast algorithm exists for the KF in 
the configuration of interest, mainly because of the structures 
of the state equation and the smoothness matrix. However, 
fast algorithm may be developed on the basis of high-order 
displacement matrices [30]. More precisely, it is easy to see 
that the displacement matrix of order 2fc + 1 (if integer) is null 
for Afc. Taking advantage of this property may result in a fast 
version of the proposed algorithm. 

However, calculation time problems are now less crucial 
than they used. The standard KS algorithm only takes 0.36 43 
to process the entire data set of Fig. [T] so, real time compu- 
tations can probably be achieved. 

'The proposed algorithm has been implemented using the computing 
environment matlab on a Personal Computer, Pentium III, with a 450 MHz 
CPU and 128 Mo of RAM. 



VI. Hyperparameters estimation 

The estimated a,„-sequence and spectra sequence depend 
on 7\/ + 4 hyperparameters: smoothness and AR orders k and 
P, power sequence r^, and two regularization parameters As 
and Ad. 

A. Power parameters 

The M parameters rj^ are needed by the proposed RegLS 
method as well as the LS and ALS procedures and the same 
empirical estimates will be used for all of them. In the cri- 
terion ([Tol l, parameters rf,, only act as weighting coefficients, 
so that the successive terms are of equivalent weight. The 
proposed empirical technique replaces the prediction error 
powers rj^ by the signal powers r„i themselves. A simple 
empirical estimate f™ = ymUm/N could be used. However, 
since the estimation variance is high for iV = 8, in practice, 
a more efficient technique consists in smoothing the sequence 
f „i. Let us note that [2] proposes a scheme which is equivalent 
in principle. 

B. Order parameters 

The proposed framework allows to estimate long AR models 
to describe various spectral shapes. Moreover, by choosing the 
maximal order P = — 1 we get rid of the difficult problem 
of model order selection. In fact, as expected and confirmed 
in Section IVII-CI as long as P is large enough it does not 
affect significantly the spectral shape. 

On the other hand, to our experience, the smoothness order 
k does not affect the spectrum sequence provided that fc ^ 0. 
So, the smoothness order is a priori tuned to fc = 1, i.e., a first 
order derivative spectra penalization. Moreover, Section [VII-CI 
also provides a quantitative sensitivity study of the spectra 
sequence w.rt. this parameter. 

C. Regularization parameters 

The problem of regularization parameter estimation within 
the proposed framework is a delicate one. It has been ex- 
tensively studied and several techniques have been proposed 
and compared [26,31-35]. The ML approach is often chosen 
within the Bayesian framework, mentioned in Remarks |2] 
and |3] The Gaussian likelihood function ([TJ and the Gaussian 
prior (|9]l together yield a Gaussian marginal law for the ob- 
served samples f{y ; Ag, Ad), i.e., the regularization parameter 
likelihood. The Hyperparameter-Co-Log-Likelihood (HCLL) 
is easily computed, for a given hyperparameter set, as a 
function of innovation vectors and covariances Rm, i.e., 
two of the KF sub-products: 



HCLL{X,,Xd) 



M 



In det R„ 



1 '■n 



ignoring constant coefficients. This expression is the gener- 
alization of a more conventional identity, available for scalar 
observations [2]. The error covariance matrix i?„j is an L x L 
matrix, L possibly ranging from L = l to L = N + P, 
according to the windowing form and model order. Since 
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i = 1 is selected in the presented computations, no specific 
algorithm has been developed for inversion nor determinant 
calculations. 

The ML estimate: 



can be computed by means of several algorithms: coordi- 
nate/gradient descent algorithm [28] or EM algorithms [36, 
37], but none of them can ensure global optimization. Here, 
the optimization stage is tackled by means of a coordinate 
descent algorithm with a golden section line search [28]. Since 
HCLL is a function of two variables only, the optimization 
stage only requires about 10 s^. 

VII. Simulation results and comparisons 

The present section assesses the effectiveness of the pro- 
posed method, compared to the usual ones by processing the 
example shown in Fig. [T] 

A. Quantitative comparison criterion 

Since the true spectrum sequence is known in the presented 
simulations, quantitative criteria are computable on the basis 
of distances between estimated spectra Sm{v) and true ones 
Sm{v), accumulated over the AI bins. Normalized distances: 

Em=l/o^ l'S'm(i^)|''dl/ 

with r = 1 and r = 2 have been computed. The normalization 
is chosen so that a null estimated spectrum results in a 100% 
error. Practically, the integrals are approximated by discrete 
summation over the frequency domain v = q/Q,q G ]Nq_i 
with Q = 1024. 



B. Tuning parameters 

1) Usual methods: Since no automatic parameters tuning 
is available for usual methods, these parameters have been 
chosen in order to produce the best distance. Moreover, we 
have checked that such a quantitative procedure finds itself in 
good agreement with the visual appreciation. 

- First of all, it is noticeable that, even for a short model, 
the non-windowed and pre-windowed methods system- 
atically yield numerous spurious peaks. The best re- 
sults have been obtained with the post-windowed forrr|| 
(double-windowed behaves similarly) so, the estimated 
spectra are of poor resolution [15, p. 225]. 

- As expected, since the true spectra show up to three 
modes, the best results have been obtained with P ~ 3 
for both LS and ALS. 

- Finally, as far as the ALS method is concerned, W — 20 
has been selected. 




Fig. 2. The Ihs and rhs figure respectively show HCLL and distance 
(L^ behaves similarly) as a function of regularization parameters (As,A(j), 
respectively read on the vertical and the horizontal axis (logjQ scaled). In 
both cases, a star (*) locates the minimum. 

2) Regularized method: The HCLL function has been com- 
puted on a fine discrete logio grid of 100 x 100 values between 
—2 and 1 for Ag and between 1 and 3 for Ad. The result is 
the HCLL sheet shown in Fig. |2l-lhs. It is fairly regular, and 
exhibits a single minimum at A"^ ~ —1.53 and A""" = 2.16. 
Moreover, Fig. |2}rhs shows the corresponding distances 
and the strikingly similar behavior of HCLL{Xs, Xd) and 
L^(As, Ad) is a strong argument in favor of the likelihood as 
a criterion for parameters tuning. 

However, it must be mentioned that a variation of on decade 
on As or Ad entails a nearly imperceptible variation in the 
estimated spectra and a fraction of percent error. This point 
is especially important for qualifying the robustness of the 
proposed method. Contrary to the choice of model order in 
the usual AR analysis, which is critical, the choice of (As, Ad) 
offers broad leeway and can be made reliably. 

Practically, the adjustment is set using the coordinate de- 
scent algorithm and Fig.|2}lhs illustrates its convergence, from 
three different starting points. 

C. Order sensitivity 

This subsection assesses the sensitivity of the method w.r.t. 
the order parameters k and P. For P 1 to P = 7 and 
for fc = to fc = 2 (step .25), we have computed the ML 
estimate (I3l1 l: 

(As"^(P, k), A^'^(P, k)) = argmini/C7iL(As, Ad, P, k) , 

and the corresponding optimal likelihood and distance 

HCLL,^,{P,k) = HCLL{XriP,k),XZiP,k),P,k) 
LSp.(P,fc) = L^iXriP,k),XTiP,k),P,k) 

They are plotted in Fig. [3] as a function of P for the several 
values of k. 

As far as the likelihood is concerned, 

• HCLLopt is a decreasing (almost linear) function of 
model order P: the ML selected order is the maximal 
one P = iV - 1 = 7. 

-A possible explanation for this rather counterintuitive fact, is that the post- 
windowed form is somewhat "self penalizing", i.e., the con'esponding criterion 
incorporates quadratic penalization terms: alnMam, where M only depends 
upon the data. 



(As"^ AD = argminifCLL(As, Ad) (31) 
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• HCLLopi does not depend on k (the four curves are 
over plotted) so that, given P the triplet (As, Ad, k) "over- 
parameterize" the likelihood and k is indifferent. 

As far as the is concerned, it still behaves similarly to the 
likelihood: it is roughly decreasing with P and not depending 
upon k. As a conclusion, the maximization of the likelihood 
w.rt. k and P does not provide any improvement and the 
recommended scheme described in Section lVI-Bl is an efficient 
one. 




Fig. 3. Optimal likeliliood HC'LLopt{P,k) (top) and distances L'^^i{P,k) 
(bottom) as a function of order P for several smoothness order k = 
0.5, 1, 1.5, and 2. 



D. Qualitative evaluation 

We have then compared the usual methods at their best 
(optimally adjusted parameters knowing the true spectra) with 
the proposed method (automatic selection of regularization 
parameters without knowledge of the true spectra). The results 
obtained by LS, ALS, and RegLS are presented in Fig. |4] A 
simple qualitative comparison with the reference Fig. [T] already 
leads to four conclusions. 

- The ML strategy provides a good value for the regular- 
ization parameters and the (and L^) distance is in 
accordance with the qualitative assessment. 

- The effect of the regularization is obvious. Estimated 
spectra are in a much greater conformity with the true 
ones. The spectrum shapes are reproduced more precisely, 
in one, two or three modes. Their positions and their 
amplitudes are correctly estimated. 

- Moreover, the spectral resolution for the ground clutter 
is strongly enhanced. It is essentially due to the coherent 
accounting for spectral and spatial continuity resulting in 
a robust non-windowed form. 

- However, it can be seen, though, that the sudden transi- 
tions at the beginning of the ground clutter is slightly 
over-smoothed. This can be expected from quadratic 
regularization and may be at least partially avoided by 
introducing non-quadratic regularization [38^0]. 

E. Quantitative evaluation 

In the non adaptive context, quantitative comparisons have 
previously been performed in [1,26]. The adaptive extension 
originally proposed by Kitagawa 8c Gersch has also been 
quantitatively assessed in [2]. 



Method 


L^ 


L^ 


Periodogram 


87.1% 


92.9% 


Best LS 


76.6% 


85.4% 


Best ALS 


66.4% 


75.5% 


ML & RegLS 


57.9% 


69.2% 



TABLE I 

Quantitative comparison of the periodogram, least squares 

METHODS and THE REGULARIZED ONE. L^ AND L? INDICATE THE 
DISTANCES BETWEEN ESTIMATED AND TRUE SPECTRA. 



For the proposed method, quantitative comparison have 
been achieved by evaluating and distances between true 
and estimated spectra. The results are listed in Table H] and 
show an L'^ improvement of about 10% form periodogram to 
best LS, 10% from best LS to best ALS and 10% from best 
ALS to the entirely automatic proposed method. 

VIII. Conclusion and perspectives 

This paper tackles short-time adaptive AR spectral esti- 
mation within the regularization framework. It proposes a 
new regularized least squares criterion accounting for spectral 
smoothness and spatial continuity. The criterion is efficiently 
optimized by a special Kalman smoother In this sense, the 
present study significantly deepens the contributions of [1,2], 
given that the latter separately address spectral smoothness and 
spatial continuity. Moreover, the proposed method is entirely 
unsupervised and it is shown that maximum likelihood regular- 
ization parameters is both formally achievable and practically 
useful. Finally, a simulated comparison study is proposed in 
the field of Doppler radars. It shows an improvement of about 
10%, comparing some usual methods at their best versus the 
entirely automatic proposed one. 

Future works will be devoted to compensate for the over- 
smoothing character of quadratic regularization in the presence 
of spatial breaks. [41] accounts for spatial continuity while 
preserving breaks by way of a non-Gaussian state model and 
extended KF algorithms. In our mind, a preferable approach 
could be to introduce non-quadratic convex penalty terms 
and to minimize the resulting criterion using descent algo- 
rithms [38,39,42]. 
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